#include "Chebyshevron.h"
#include "FourierTron.h"
#include "Polynomicon.h"
#include "Mobius.h"
#include <iostream>


using namespace std;
using namespace SeamStress;


Chebyshevron::Chebyshevron(Mobius* mob, vector<Seamstress*>* ss, unsigned long int nt, mpfr_prec_t p) : pins(this, ss)
{
  if(nt==0){nthreads=1;}
  else{nthreads=nt;}
  if(nthreads > ss->size()){nthreads = ss->size();}
  
  prec = p;
  
  mobius = mob;
  
  tempval1.set(0,0,prec);
  tempval2.set(0,0,prec);
  one.set(1,0,prec);
  negone.set(-1,0,prec);
  
  threadvec1.assign(nthreads, tempval1);
}


Chebyshevron::~Chebyshevron()
{
  
}


void Chebyshevron::valsFromCoefs(const std::vector<cxfloat>& coefs, std::vector<cxfloat>& vals)
{
  unsigned long int thread_count=nthreads;
  if(coefs.size()<thread_count){thread_count=coefs.size();}
  
  tempval1.set(0,0);
  tempvec1.resize(coefs.size(),tempval1);
  for(unsigned long int i=0;i<coefs.size();i++){tempvec1[i].set(coefs[i]);}
  
  mobius->transform(one, one, negone, one, tempvec1, tempval1);
  
  tempval1.set(0,0);
  unsigned long int N=tempvec1.size();
  tempvec1.resize((N<<1)-1, tempval1);
  for(unsigned long int i=N;--i;)
  {
    __builtin_prefetch(&(tempvec1[i-1]), 0, 3);
    __builtin_prefetch(&(tempvec1[(i-1)<<1]), 1, 3);
    tempvec1[i<<1].set(tempvec1[i]);
    if((i&1)==1){tempvec1[i].set(tempval1);}
  }
  mobius->transform(negone, one, one, one, tempvec1, tempval1);
  
  long int scale = -((N-1)<<1);
  unsigned long int M = tempvec1.size() - 1;
  if((N&1)==1)
  {
    for(unsigned long int i=0;i<M;i++)
    {
      __builtin_prefetch(&(tempvec1[i+1]), 1, 3);
      tempvec1[i].mul_power_2(scale);
    }
    tempvec1.back().mul_power_2(scale);
  }
  else
  {
    for(unsigned long int i=0;i<M;i++)
    {
      __builtin_prefetch(&(tempvec1[i+1]), 1, 3);
      tempvec1[i].mul_power_2(scale);
      tempvec1[i].neg();
    }
    tempvec1.back().mul_power_2(scale);
    tempvec1.back().neg();
  }
  
  tempval1.set(0,0);
  M = 1;
  thread_var2 = 0;
  N<<=1;
  while(M<N){M<<=1;thread_var2++;}
  tempvec1.resize(M, tempval1);
  
  mobius->getPolynomicon()->getFourierTron()->calculateTwiddle(N);
  
  //t = w*exp(i*pi/M)
  //w = t*exp(-i*pi/M)
  //w^k = t^k*exp(-i*pi*k/M)
  //M = 2^thread_var2
  
  thread_tot = thread_count;
  thread_var1 = M;
  pins.sewStraight(&Chebyshevron::scaleVarByUnityRoot, thread_tot);
  
  
  
  mobius->getPolynomicon()->getFourierTron()->CT_FFT(tempvec1, M);
  
  N = (M>>1);
  vals.resize(N, tempval1);
  
  thread_var1 = N-1;
  thread_tot = thread_count;
  tempvecpointer1 = &vals;
  thread_var3 = coefs.size() - 1;
  pins.sewStraight(&Chebyshevron::finalizeVals, thread_tot);
}


void Chebyshevron::scaleVarByUnityRoot(void* arg)
{
  unsigned long int w = (*((unsigned long int *)arg));
  
  unsigned long int i = w;
  unsigned long int max = thread_var1 - thread_tot;
  
  for(i=w;i<max;i+=thread_tot)
  {
    __builtin_prefetch(&(tempvec1[i+thread_tot]), 1, 3);
    mobius->getPolynomicon()->getFourierTron()->twiddleVal(thread_var2, i, threadvec1[w]);
    tempvec1[i].mul(threadvec1[w]);
  }
  mobius->getPolynomicon()->getFourierTron()->twiddleVal(thread_var2, i, threadvec1[w]);
  tempvec1[i].mul(threadvec1[w]);
}


void Chebyshevron::finalizeVals(void* arg)
{
  unsigned long int w = (*((unsigned long int *)arg));
  
  unsigned long int max = thread_var1 - thread_tot;
  unsigned long int i=w;
  unsigned long int k=0;
  unsigned long int N = ((tempvec1.size())<<1);
  
  for(i=w;i<=max;i+=thread_tot)
  {
    __builtin_prefetch(&(tempvec1[i+thread_tot]), 0, 3);
    __builtin_prefetch(&((*tempvecpointer1)[i+thread_tot]), 1, 3);
    k = N - ((((i<<1) + 1)*thread_var3))%N;
    mobius->getPolynomicon()->getFourierTron()->twiddleVal(thread_var2, k, (*tempvecpointer1)[thread_var1-i]);
    (*tempvecpointer1)[thread_var1-i].mul(tempvec1[i]);
  }
  k = N - ((((i<<1) + 1)*thread_var3))%N;
  mobius->getPolynomicon()->getFourierTron()->twiddleVal(thread_var2, k, (*tempvecpointer1)[thread_var1-i]);
  (*tempvecpointer1)[thread_var1-i].mul(tempvec1[i]);
}


void Chebyshevron::coefsFromVals(const std::vector<cxfloat>& vals, std::vector<cxfloat>& coefs)
{
  //check that vals.size() is a power of two
  try{if((vals.size() & (vals.size()-1)) != 0){throw (char*)("size of input must be a power of two");}}
  catch(char* str){cout<<"Exception from Chebyshevron::coefsFromVals : "<<str<<endl;throw;}
  
  coefs.resize(vals.size(), tempval1);
  
  unsigned long int thread_count=nthreads;
  if(vals.size()<(1+thread_count)){thread_count=vals.size()-1;}
  
  tempvec1.resize((vals.size())<<1, tempval1);
  for(unsigned long int i=0;i<vals.size();i++)
  {
    tempvec1[i].set(vals[vals.size() - 1 - i]);
    tempvec1[tempvec1.size() - 1 - i].set(tempvec1[i]);
  }
  thread_var1 = vals.size() - 1;
  thread_tot = thread_count;
  thread_var2 = 1;
  unsigned long int M = 1;
  while(M<vals.size()){M<<=1;thread_var2++;}
  mobius->getPolynomicon()->getFourierTron()->calculateTwiddle(M);
  pins.sewStraight(&Chebyshevron::initializeVals, thread_tot);
  mobius->getPolynomicon()->getFourierTron()->inverseCT_FFT(tempvec1, tempvec1.size());
  
  thread_var1 = tempvec1.size();
  thread_var3 = thread_var1<<1;
  pins.sewStraight(&Chebyshevron::scaleVarByUnityRootInverse, thread_tot);
  unsigned long int N = vals.size();
  tempvec1.resize(((N - 1)<<1) + 1);
  M = tempvec1.size() - 1;
  long int scale = (N - 1)<<1;
  if((N&1)==1)
  {
    for(unsigned long int i=0;i<M;i++)
    {
      __builtin_prefetch(&(tempvec1[i+1]), 1, 3);
      tempvec1[i].mul_power_2(scale);
    }
    tempvec1.back().mul_power_2(scale);
  }
  else
  {
    for(unsigned long int i=0;i<M;i++)
    {
      __builtin_prefetch(&(tempvec1[i+1]), 1, 3);
      tempvec1[i].mul_power_2(scale);
      tempvec1[i].neg();
    }
    tempvec1.back().mul_power_2(scale);
    tempvec1.back().neg();
  }
  
  mobius->transform(negone, one, one, one, tempvec1, tempval1);
  scale = -scale;
  for(unsigned long int i=0;i<M;i++)
  {
    __builtin_prefetch(&(tempvec1[i+1]), 1, 3);
    tempvec1[i].mul_power_2(scale);
  }
  tempvec1.back().mul_power_2(scale);
  
  N = (tempvec1.size() - 1)>>1;
  for(unsigned long int i=1;i<N;i++)
  {
    __builtin_prefetch(&(tempvec1[i+1]), 1, 3);
    __builtin_prefetch(&(tempvec1[(i+1)<<1]), 1, 3);
    tempvec1[i].set(tempvec1[i<<1]);
  }
  tempvec1[N].set(tempvec1[N<<1]);
  tempvec1.resize(N+1);
  
  mobius->transform(one, negone, one, one, tempvec1, tempval1);
  for(unsigned long int i=0;i<coefs.size();i++){coefs[i].set(tempvec1[i]);}
  N = vals.size();
  M = coefs.size() - 1;
  scale>>=1;
  if((N&1)==1)
  {
    for(unsigned long int i=0;i<M;i++)
    {
      __builtin_prefetch(&(coefs[i+1]), 1, 3);
      coefs[i].mul_power_2(scale);
    }
    coefs.back().mul_power_2(scale);
  }
  else
  {
    for(unsigned long int i=0;i<M;i++)
    {
      __builtin_prefetch(&(coefs[i+1]), 1, 3);
      coefs[i].mul_power_2(scale);
      coefs[i].neg();
    }
    coefs.back().mul_power_2(scale);
    coefs.back().neg();
  }
}


void Chebyshevron::initializeVals(void* arg)
{
  unsigned long int w = (*((unsigned long int *)arg));
  
  unsigned long int max = thread_var1 - thread_tot;
  unsigned long int i=w;
  unsigned long int k=0;
  unsigned long int N = ((tempvec1.size())<<1);
  
  for(i=w;i<=max;i+=thread_tot)
  {
    __builtin_prefetch(&(tempvec1[(i+thread_tot)]), 1, 3);
    __builtin_prefetch(&(tempvec1[tempvec1.size() - (i+thread_tot)]), 1, 3);
    
    k = ((((i<<1) + 1)*thread_var1))%N;
    mobius->getPolynomicon()->getFourierTron()->twiddleVal(thread_var2, k, threadvec1[w]);
    tempvec1[i].mul(threadvec1[w]);
    k = N - k;
    mobius->getPolynomicon()->getFourierTron()->twiddleVal(thread_var2, k, threadvec1[w]);
    tempvec1[tempvec1.size() - 1 - i].mul(threadvec1[w]);
  }
  k = ((((i<<1) + 1)*thread_var1))%N;
  mobius->getPolynomicon()->getFourierTron()->twiddleVal(thread_var2, k, threadvec1[w]);
  tempvec1[i].mul(threadvec1[w]);
  
  k = N - k;
  mobius->getPolynomicon()->getFourierTron()->twiddleVal(thread_var2, k, threadvec1[w]);
  tempvec1[tempvec1.size() - 1 - i].mul(threadvec1[w]);
}


void Chebyshevron::scaleVarByUnityRootInverse(void* arg)
{
  unsigned long int w = (*((unsigned long int *)arg));
  
  unsigned long int i = w;
  unsigned long int max = thread_var1 - thread_tot;
  
  for(i=w;i<max;i+=thread_tot)
  {
    __builtin_prefetch(&(tempvec1[i+thread_tot]), 1, 3);
    mobius->getPolynomicon()->getFourierTron()->twiddleVal(thread_var2, thread_var3 - i, threadvec1[w]);
    tempvec1[i].mul(threadvec1[w]);
  }
  mobius->getPolynomicon()->getFourierTron()->twiddleVal(thread_var2, thread_var3 - i, threadvec1[w]);
  tempvec1[i].mul(threadvec1[w]);
}


void Chebyshevron::chebyPoly(unsigned long int N, refloat_vec& coefs)
{
  //we want to interpolate between a power of 2 # of points, even if N is not a power of 2
  thread_var3 = N-1;
  thread_var3 = (thread_var3 >> 1) | thread_var3;
  thread_var3 = (thread_var3 >> 2) | thread_var3;
  thread_var3 = (thread_var3 >> 4) | thread_var3;
  thread_var3 = (thread_var3 >> 8) | thread_var3;
  thread_var3 = (thread_var3 >> 16) | thread_var3;
  thread_var3 = (thread_var3 >> 32) | thread_var3;
  thread_var3++;
  tempvec2.resize(thread_var3, tempval1);
  //   thread_var2 = log_2(N)
  unsigned long int M = thread_var3>>1;
  thread_var2 = 0;
  while(M)
  {
    thread_var2++;
    M >>= 1;
  }
  thread_var1 = N;
  
  thread_tot=nthreads;
  if(N<thread_tot){thread_tot=N;}
  mobius->getPolynomicon()->getFourierTron()->calculateTwiddle(N<<1);
  pins.sewStraight(&Chebyshevron::chebyPolyInit, thread_tot);
  
  coefsFromVals(tempvec2, tempvec3);
  refloat reval(0,prec);
  coefs.resize(N, reval);
  for(unsigned long int i=0;i<N;i++){tempvec3[i].get_real(coefs[i]);}
  if((N&1) == 0)
  {
    for(unsigned long int i=0;i<N;i++){coefs[i].neg();}
    for(unsigned long int i=0;i<N;i+=2){coefs[i].set(0);}
  }
  else
  {
    for(unsigned long int i=1;i<N;i+=2){coefs[i].set(0);}
  }
}


void Chebyshevron::chebyPolyInit(void* arg)
{
  unsigned long int w = (*((unsigned long int *)arg));
  
  //n = N
  //x_k = cos( pi/2 * (2*k+1)/n ), k in [0, N)
  //T_(N-1)(x_k) = cos( (N-1) * pi/2 * (2*k+1)/n )
  
  
  for(unsigned long int i=w;i<thread_var3;i+=thread_tot)
  {
//     unsigned long int k = thread_var1 - i;
    unsigned long int k = i;
    unsigned long int j = (2*k + 1)*(thread_var1 - 1);
    j = j%(thread_var3*4);
    j = thread_var3*4 - j;
    
    mobius->getPolynomicon()->getFourierTron()->twiddleVal(thread_var2+1, j, tempvec2[i]);
    tempvec2[i].set_imag(0);
  }
}



